/*  sifbdfct.cpp  by K.Tsuru  */
// function ID = 4106   BRADIX
/********************************************************************
SInteger class
It provides the factorial n! using the partial products.
It is a few times faster than BFact(), because the FFT multiplication
can be used.
*********************************************************************/
#ifndef SN_H
#include "sn.h"
#endif

/*********************************************************************
sub-function
It provides a product of continuous integers begining with f(<=m).
r = f(f+1)....g
If the processing ended(g==m) it returns zero, else "g+1" i.e. the
next integer.
**********************************************************************/
static ulong MultContNum(ulong f, ulong m, uint fig, SInteger& r){
    r.SetLong(f);
    if(f == m) return 0;

    ulong i = f+1;
    int end = 0;
    while((r.Head() < fig) && !end){
        IsMult(r, i, r);
        if(i == m) end = 1;
        i++;
    }
    return  end ? 0 : i;
}
// SNBlock version n!
SInteger BdFact(ulong n){
    long fig = long( (double)Stirling(n)/log10((double)BRADIX) ) +1;//the number of figures in BRADIX ver. 2.17
    if(fig >= (long)SNManager::SNMaxSize(SNManager::BIN_INT)){
        SNManager::SetError(SNManager::OUT_OF_RANGE, "BdFact", 4106);
    }

    SInteger m; // m = 1
    if(n < 2LU){    m.SetLong(1L);  return m;   }

    uint size_mul = 2u;//The speed depends on this value.
    uint size = size_mul*m.FFTMinSize(), s2 = uint(fig/(long)(size-2u));
    s2 = ceilpow2(s2+1u);   //+1u considering s2 == 0
    SNBlock <SInteger> w(s2);   // work area
    ulong i =2LU;

// step 1 : It makes partial products whose figure is "size" or lower.
    uint s, k = 0;
    while(i && (k < s2)){
        w[k].FigureAlloc(size, 0);
        i = MultContNum(i, n, size-2u, w[k++]);
    }
    if( (i == 0) && (k == 1u) ) return w[0];
    for(;k < s2; k++) w[k].SetSmall(1);//It fills the residual part of w[] with one.
// step 2 : It composes the partial products.
    s = s2;
    while(s >= 2u){
        for(k = 0; k <= s-2u; k += 2u){
        //When w[k] == 1, w[k+1] == 1.
            if(w[k].IsOne() && w[k+1u].IsOne()) w[k/2u].SetSmall(1);
            else if(w[k+1].IsOne()) w[k/2u] = w[k];
            else                    w[k/2u] = w[k+1u]*w[k]; //The FFT multiplication is used.
        }
        s /= 2u;
        //It frees the memory.
        for(k = s; k < 2u*s; k++) w[k].SizeZero();
    }
    while(i){   //recidual part exists? perhaps not
        i = MultContNum(i, n, size-2u, m);
        w[0] *= m;
    }
    return w[0];    // SNBlock class's destructor frees the memory.
}